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Abstract 

We relate progress in statistical mechanics, both at and far from equilibrium, to advances in the 
theory of dynamical systems. We consider computer simulations of time-reversible deterministic 
chaos in small systems with three- and four-dimensional phase spaces. These models provide us with 
a basis for understanding equilibration and thermodynamic irreversibility in terms of Lyapunov 
instability, fractal distributions, and thermal constraints. 
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I. INTRODUCTION 


Nonequilibrium Molecular Dynamics and Dynamical Systems Theory have been our main 
research interests for about 50 years, the same period over which Moore’s Law has described 
the growth of our primary tool, computation. In 1959 thermodynamic information was 
mainly gleaned from series expansions of pressure in powers of the density and integral 
equations for the pair distribution function. That was the year when Berni Alder and Tom 
Wainwright described a new simulation method- now called “molecular dynamics” in their 
prescient Scientific American article “Molecular Motions” : 

“One of the aims of molecular physics is to account for the bulk properties of 
matter [ pressure P , temperature T, energy E, ... ] in terms of the behavior of 
its particles. High-speed computers are helping physicists realize this goal.” 

At that time simulating the motion of a few hundred particles presented a computational 
challenge. Today’s biomolecule simulations model at least many thousand and perhaps a few 
million atomistic degrees of freedom. After several Nobel prizes^ this molecular dynamics 
method is familiar textbook material while the virial series for the pressure and the pair- 
distribution integral equations keep company with the dinosaurs. 

During this same period our understanding of dynamical systems ( flows described by a 
few nonlinear ordinary differential equations ) has undergone explosive growth. Ed Lorenz’ 
three-equation Butterfly Attractor is a clearcut demonstration of “chaos”, the exponen¬ 
tial “Lyapunov instability” often found in systems of three or more ordinary differential 
equations. The Lyapunov spectrum of exponential growth and decay rates provides a topo¬ 
logical description of evolving phase-space densities. The discovery that time-reversible flow 
equations can describe irreversibility through the formation of fractal strange attractors fur¬ 
nished a new geometric interpretation of the Second Law of Thermodynamics in terms of 
an underlying reversible mechanics^. 

The correspondence between manybody molecular dynamics and the concepts developed 
in dynamical systems theory involves five key ideas : 

[ 1 ] Simulating nonequilibrium systems requires a new nonequilibrium molecular dynamics 
which, unlike Hamiltonian mechanics, includes thermodynamic control variables. 

[ 2 ] These control variables, such as thermostats or ergostats, can provide ergodic equi¬ 
librium dynamics, replicating Gibbs’ canonical distribution. 
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[ 3 ] Away from equilibrium this same approach, while time-reversible, can promote and 
maintain noneqnilibrium steady states. 

[ 4 ] Despite the time-reversible nature of the nonequilibrium flow equations the resulting 
phase-space description is dissipative, on average, and generates multifractal attractors. 

[ 5 ] The multifractal nature of nonequilibrium steady states confirms their rarity and 
provides a mechanical explanation of the Second Law of Thermodynamics. 

The dynamical systems approach to irreversible processes continues to provide new insight 
into both equilibrium and far-from-equilibrium flows^. Our intent here is to illustrate this 
insight by the exploration of the simplest possible dynamical models for nonequilibrium 
steady states. We begin with the Galton Board problem^ -7 , a steady field-driven flow with 
impulsive hard-disk collisions. We continue, and then conclude, with a variety of generalized 
harmonic oscillator problems^. These illustrate heat flow and ergodic fractal formation with 
just three ordinary differential equations. The hard-disk Galton Board and the generalized 
conducting-oscillator problems display all of the key ideas linking manybody mechanics to 
small-system analogs. 

The plan of this work is as follows. We first review the Galton Board problem and use 
that example to illustrate the fractal attractors generated by time-reversible nonequilibrium 
steady states. We visualize these attractors through two-dimensional cross sections of their 
three-dimensional phase-space distributions. The Galton Board is one of the simplest chaotic 
problems. It is deterministic and ergodic in its three-dimensional phase space. The ergodicity 
is enabled by the ( exponential ) Lyapunov instability of its hard-disk collisions. 

We then explore ergodicity ( dynamical access to all phase-space states ) for smoothly- 
continuous harmonic-oscillator problems, at and away from thermal equilibrium. Many of 
the nonequilibrium versions of oscillator problems provide dissipative strange attractors in 
just three or four phase-space dimensions. We point out some useful numerical techniques 
for exploring the boundary between chaos and regularity and discuss the possibility of nu¬ 
merical implementations of Liouvillc’s phase-space flow equations. Finally we tie together 
these simple microscopic example problems to their real-world analogs in macroscopic ther¬ 
modynamics and computational fluid mechanics. 
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FIG. 1: Description of collisions in the Galton Board in terms of the angles a and (5 . At a collision 
the radial velocity, — cos(/3) changes sign while the tangential velocity sin(/3) is unchanged. 

II. THE GALTON BOARD ERGODIC, TIME-REVERSIBLE, DISSIPATIVE 

Our goal throughout is to connect dynamics, statistical mechanics, and nonequilibrium 
systems with the simplest possible examples. The Galton Board models Sir Francis Galton’s 
lecture table probability demonstration based on the chaotic motion of particles introduced 
at the top and in the center of a fixed lattice of scatterers. In our idealized mechanical 
steady-state model this field-driven motion occurs at constant speed. To implement this 
idea the field’s acceleration is moderated by a deterministic time-reversible “thermostat” 
force acting parallel to the particle’s velocity. The resulting trajectory is isokinetic and 
continuously dissipates field energy as heat. The field is the source of energy. The heat 
reservoir represented by a thermostat force is the compensating heat sink. The overall 
description of this isothermal ( constant kinetic energy ) Galton Board model system is the 
prototypical simplest mechanical example of a time-reversible nonequilibrium steady state. 
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FIG. 2: Field Dependence of Galton Board Collisions. Each collision corresponds to a single 
point [ 0 < a < tt with — 1 < sin(/3) < +1 ] so that each of these four sample phase-space cross 
sections illustrates 300,000 successive collisions. At zero field strength the distribution of points is 
completely uniform with a constant density of points. 

The computation of its dynamics involves solving four coupled ordinary differential equations 
for the (x,y) location of the falling particle and its velocity ( p x ,P y ) and is punctuated by 
hard-disk scatterer collisions. For simplicity the falling particle has unit mass and speed. 
We choose the accelerating field parallel to the y axis, which is perpendicular to one set of 
rows of scatterer particles, as shown in Figure 1. 

The resulting diffusive motion through such a periodic array of scatterers is easy to 
program, particularly if the scatterers are motionless “hard disks”. One simply integrates 
the motion equations for ( x , y,p x ,P y ) until the moving particle finds itself “inside” a scatterer 
( | r | < 1/2 ) . Then the dynamics is returned to the previous ( x,y ) coordinate set. There 
the sign of the radial velocity is changed from negative to positive, and the integration 
is continued. We have found that Runge-Kutta integration is the simplest useful method 
for generating trajectories between collisions. The alternative analytic approach^, though 
feasible, is cumbersome. 
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In between collisions, the equations of motion, with the field —E in the y direction, are : 


x = p x ; y = p y ; p x = -C Vx ; i> y = -E - (p y ■ 

We choose the linear “frictional force” —(p, to enforce an isokinetic constant-speed con¬ 
straint. This linear force is the simplest choice. It also follows directly from Gauss’ Principle 
of Least Constraint^, where the constraint condition controls the kinetic energy, for which 
we choose 

K = (p 2 , + pl )/2 = (1/2); K = 0. 

The linear constraint force, — (p , is sufficient to satisfy the isokinetic condition : 

K = 0 = P X j) X + PyPy = -(Pi - Ep y - Qjy - > C = ~ E Py ' 

A unit cell, within which the motion occurs, is pictured in Figure 1 . It is convenient 
to think of the moving particle as a mass point and the scatterers as fixed particles of unit 
diameter. At each collision ( defined by the two angles a and [3 shown in Figure 1 ), the 
radial component of the velocity [ — cos (/3) } is reversed from negative to positive and the 
motion is continued with the resulting post-collision values of { p x ,p y } . By choosing a 
scatterer density of four fifths the maximum close-packed density we avoid the possibility of 
a ballistic collisionless trajectory. The inevitable scatterer collisions make possible ( and for 
moderate field strengths, inevitable ) diffusive piecewise-continuous trajectories punctuated 
by a series of scatterer collisions. 

Three distinct types of solution result, conservative, dissipative, and periodic, with the 
type determined by the initial condition and the field strength E y . For zero held the mo¬ 
tion is ergodic and conservative, obeying the equilibrium version of Liouville’s Theorem, 
f(x,y,p x ,p y ) = 0 . That is, all conceivable collision types do occur, and with a uniform 
probability when plotted in the [ a, sin(/3) ] plane. As the held is increased it becomes ap¬ 
parent that the distribution of collisions, though becoming nonuniform, remains “ergodic”, 
with nonzero probability everywhere. We use the word “ergodic” as equivalent to the Ehren- 
fests’ “quasiergodic” notion of a dynamics that eventually comes arbitrarily close to every 
point in the distribution. At relatively high values of the accelerating held things can be 
different. Trajectories can become trapped in stable periodic orbits, some conservative and 
some dissipative so that the motion is no longer ergodic. See, for an obvious example, the 
“hole” in the distribution for E = 4 shown in Figure 2 . 
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Quantitative investigation of the two-dimensional phase-space cross sections illustrated 
in Figure 2 reveals that the densities of points in the vicinity of ( 1 ) a randomly chosen 
point oc r Dl and of ( 2 ) a collisional point oc r D2 are different, with 2 > D\ > D 2 > 1 
. These power laws define^ ( 1 ) the “information dimension” and ( 2 ) the “correlation 
dimension” of the various fractal cross sections. Both of these two fractional dimensions 
vary with field strength. 

Despite this singular fractal behavior there is no problem computing the mean vertical 
current and the equivalent conductivity. New phenomena appear at field strengths somewhat 
higher than E = 3 : stable sequences of repeated collisions begin to occur. In the full three- 
dimensional phase space, which includes an additional time dimension for the free-flight 
portion of the trajectories, the stable sequences are described by regular conservative tori 
in [ a,sin(/3),i ] space or by dissipative limit cycles in which field energy is absorbed by the 
time-reversible friction force — (p . Examples of both types are shown in Figure 3. The 
family of two-bounce horizontal orbits shown at the left has no net current. In contrast, 
the ten-bounce orbit on the right is strongly dissipative with a net downward current. This 
motion generates a one-dimensional limit-cycle orbit in the three-dimensional phase space 
and would be represented by five zero-dimensional dots in a collisional cross-section picture 
of the type shown in Figure 2 . 

Let us summarize our findings from this simple nonequilibrium steady-state problem. 
The Galton Board is deterministic, time-reversible, and dissipative. With the field “off” the 
motion is ergodic - it comes arbitrarily close to any collision type from headon [ sin(/3) = 0 ] 
to glancing [ sin(/3) = ±7r ] and anywhere from the top ( a — it ) to the bottom ( a = 0 ) of 
a scatterer. This ergodicity provides a direct connection between Newton’s dynamical and 
Gibbs’ statistical treatments of mechanical systems. 

With the field “on” the dynamics becomes fractal, though still ergodic for moderate field 
strengths. The phase-space description becomes a fractional-dimensional representation of 
collisions which is singular, nonuniform, and dissipative. The dissipation reflects the conver¬ 
sion of gravitational work into heat through the mechanism of a heat reservoir. Although 
the motion and the motion equations are perfectly time-reversible, the typical field-driven 
case is at the same time dissipative. The conversion of field energy Ey to extracted heat is 
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FIG. 3: Conservative Tori at E = 4 and a Dissipative Limit Cycle at E = 6 . The two-bounce 
trajectories shown at the left occur in the prominent “hole” seen in the E = 4 collision plot of 
Figure 2 . The periodic orbit at the right corresponds to five zero-area points in the [ a, sin(/3) ] 
representation of Figure 2 . 

imposed by the friction coefficient ( : 

(( p *) = (2(K) = (() = (S/k)>0. 

Dividing the dissipated heat (p 2 by the temperature, T = ( p 2 /k ) shows that the friction 
coefficient is also equal to the instantaneous irreversible entropy production, ( = ( S/k ) 
where k is Boltzmann’s constant. In our numerical work we set it equal to unity, k = 1 . 

Because the Galton Board distributions are fractal, with zero-area cross sections [ having 
fractal dimensionalities less than two ] the random-sampling probability of finding a point 




[ a, sin(/3) ] on the strange attractor is precisely zero. If one attempts to define a limiting 
probability density in the cross sections by counting points in small cells, he soon discovers 
that the density does not have a small-cell limiting value. Instead it diverges as a fractional 
power of the cell size. This finding was both surprising and illuminating in 1987-^. Since then 
it has turned out that such fractal attractors are typical representations of nonequilibrium 
steady states, and not just for small systems. Manybody simulations of time-reversible shear 
flows and heat flows likewise provide strange attractors with fractional dimensionalities less 
than that of the phase spaces in which they are embedded—. These attractors are “strange” 

[ fractal ] and “chaotic” [ because small perturbations on them grow exponentially fast ], 
despite the continuous equations that generate them and despite their zero-volume attractive 
nature. The dimensionality loss exhibited by attractors increases in an irregular manner with 
the field-induced departure from equilibrium. 

These fractal distributions are fully consistent with the Second Law of Thermodynamics. 
That Law declares that only the dissipative forward-in-time versions of the nonequilibrium 
trajectories are observable. The time-reversed versions of the dynamics - unphysical tra¬ 
jectories which convert heat-reservoir energy into work - are both mechanically unstable 
and computationally unobservable. Nonequilibrium dissipative trajectories seek out fractal 
attractors when followed forward in time. Reversed trajectories make up an unstable unob¬ 
servable phase-space repellor, a fractal phase-space object which repels rather than attracts 
nearby trajectories and is unstable to perturbations. Picturing such a repellor in [ a, sin(/3) ] 
space is easy. Simply reflect the fractal attractor objects of Figure 2 about their horizontal 
center lines. This changes the sign of the velocity at each collision and is equivalent to 
picturing the motion backward in time. 

The ubiquitous fractal nature of nonequilibrium steady states, singular everywhere, in¬ 
dicates the difficulty inherent in attempting their mathematical description. Although the 
motion of the hard-disk-scatterer Galton Board problem is ergodic for moderate fields, with 
all states accessible, typically, for smoothly-continuous potentials there are nonergodic situ¬ 
ations for simple mechanical systems. The mechanical treatment of theoretical models with 
smooth potential minima is complicated by the complex phase-space structure of Hamil¬ 
tonian chaos - “islands” ( the two-dimensional cross sections of three-dimensional tori ), 
chains of islands, . . . , and the endless details of structure on all scales, from large to the 
microscopic and to the unobservably smalLA 
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This unsettling complex situation can be avoided, or simplified, by the judicious use of 
thermostating forces. The best problem area for such explorations is the one-dimensional 
harmonic oscillator, the prototype of smooth Hamiltonian systems. We will see how a ther- 
mostated oscillator can be modeled so as to avoid the infinitely-many nonergodic solutions 
of Hamiltonian mechanics while providing insight into irreversible processes described by 
simple phenomenological laws. Let us turn next to the thermostating of that simplest case, 
a single harmonic oscillator with a specified temperature T rather than a fixed Hamiltonian 
energy E . 

III. NOSE SEEKS GIBBS’ CANONICAL ENSEMBLE THROUGH CHAOS 

Willard Gibbs invented his “canonical” ( in the sense of “simplest” or “prototypical” ) 
ensemble in order to link microscopic phase-space dynamics to macroscopic temperature and 
thermodynamics. His canonical ensemble collects together all the energy states accessible 
to a system in contact with a heat reservoir at a temperature T. The relative weight of 
each such state in the ensemble is the familiar Maxwell-Boltzmann weighting proportional 
to e~ E / kT . Here k is Boltzmann’s constant. 

For simplicity we focus on the application of Gibbs’ ensembles to the one-dimensional 
harmonic oscillator. With the mass and force constant and Boltzmann’s constant all set 
equal to unity Gibbs’ canonical weighting of the oscillator states is the familiar Gaussian 
distribution, a probability density for q and p : 

f(q,P ) = (2irT)-‘ e -^e-^ . 

A common textbook rationalization of the canonical distribution is to imagine that the 
members of an ensemble of systems are weakly coupled to one another. The coupling 
permits energy to be exchanged among the systems, resulting in Gibbs’ maximum-entropy 
canonical distribution. Shuichi Nose developed a much simpler picture in which a single 
system is coupled dynamically to a computational heat reservoir in such a way that a long¬ 
time average of that system’s properties is identical to the canonical average. Let us describe 
this idea in the context of the one-dimensional harmonic oscillator. 
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A. Nose’s Canonical Mechanics 


By 1984 Shuichi Nose had documented his modification of Hamilton’s constant-energy 
dynamics in two ground-breaking paper&hy^. His new dynamics was formally consistent 
with Gibbs’ constant-temperature canonical distribution. For the oscillator problem at the 
temperature T , the simplest form of Nose’s novel Hamiltonian, now with two degrees of 
freedom, # = 2 , rather than one, has the form : 

2n NoS e(q,P,s,C) = (p/s) 2 + q 2 + #Tln(.s 2 ) + ( 2 . 

The added thermostat degree of freedom, s and its conjugate momentum £ = p s , along with 
the usual (q,p) pair describes the canonical oscillator problem with four ordinary differential 
equations rather than the usual two : 

q = (p/s 2 ) ; p = - q - s = C 5 C = ( P 2 /s 3 ) - (#T/s) . [ Nose ] 

Nose carried out the straightforward and tedious algebra necessary to show that this 
approach can be made consistent with Gibbs’ canonical distribution. Three steps were 
involved in his demonstration : 

[ 1 ] “Time Scaling”: ( d/dt ) = ' —» s(d/dt) = s ' ; 

[ 2 ] redefine momentum: (p/s) —> p ; 

[ 3 ] redefine degrees of freedom: —» ff — 1 . 

More than a decade passed before Carl Dettmann simplified this approach- 4,1 —. He showed 
that starting out with a scaled Hamiltonian and setting it equal to zero , 

Ffo = sF^Nose — d j 

produces a dynamics identical to Nose’s three-step approach without the need for an explicit 
time-scaling step. 

B. Nose-Hoover Canonical Mechanics 

The even simpler “Nose-Hoover” version of Nose’s approach™ eliminates all three of these 
steps as well as the extraneous variable s . It is based on the application of Liouville’s phase- 
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FIG. 4: Trajectory intersections with the q = 0, p = 0, and £ = 0 planes are shown. The points 
in these cross sections all correspond to penetrations by a single chaotic trajectory with the initial 
conditions (q,p,() = (0,5,0) . The temperature is unity for these three equilibrium Nose-Hoover 
oscillator cross sections. 

space continuity equation to the oscillator equations of motion augmented by the definition 
of a time-dependent friction coefficient £ : 

q = p;p = ~q-Cp; £ = [ (p 2 /T) - 1 ] —> ( p 2 ) = T . [ NH ] 

This friction coefficient acts as a “thermostat”, steering the instantaneous temperature p 2 
toward the target thermostat temperature T . We can verify that this three-equation model 
is consistent with the canonical distribution for q and p augmented by a Gaussian distribution 
for £ : 

f(q,p,C) X (2vrpT = e ~i 2 /^ e - P y2T e -^/2 
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To show this we evaluate the time-rate-of-change of the probability density f(q,p, () at a 
fixed location ( q, p, C ) hi the three-dimensional phase space where the local flow velocity is 

v = ( q,p, C ) : 

(df/dt) = -V • (fv) = -f (dp/dp) - ( df/dq)q - (df/dp)p - (df/d()( = 

fC ± (■ qf/T)(p ) + (pf/T)(-q - C p) + (C/)[ (p 2 /^) «- 1 ] = o . 

Because (df/dt) vanishes everywhere the Nose-Hoover equations are consistent with Gibbs’ 
canonical distribution. 

Although the smooth and simple three-dimensional Gaussian distribution is the exact 
stationary solution of the Nose-Hoover motion equations, the new dynamics conceals an 
intricate complexity connected with “chaos”, the exponential sensitivity of calculated 
trajectories to perturbations of the initial conditions. The three [ NH ] motion equations 
included here are just the necessary minimum for chaos. But these necessary three are 
not necessarily sufficient, as is easily revealed by a numerical exploration of the oscillator’s 
phase space distribution f(q,p, () . Whether or not there is chaos depends, in a highly- 
complicated way, upon the initial conditions. 

To outline the chaos’ profile, let us advance a trajectory beginning at a known chaotic 
initial condition such as ( q,p, f ) = ( 0, 5, 0 ) or ( 3, 3, 0 ) . By plotting any two of the 
three variables { q,p, C } just as the third passes through zero, three separate cross sections 
of the chaotic sea are revealed. Figure 4 shows the sequence of 768,460 such successive 
crossing points following from the initial condition { q,p, ( } = { 0.00, 5.00, 0.00 } using one 
billion fourth-order Runge-Kutta timesteps with dt = 0.001 . We see from Figure 4 that 
only about six percent of the three-dimensional stationary Gaussian measure makes up the 
chaotic sea. Any trajectory started in the sea cannot leave and eventually explores all of it. 

The remaining phase space is occupied privately by concentric tori enclosing stable pe¬ 
riodic orbits. The simplest such orbit for the Nose-Hoover oscillator is shown in Figure 5. 
The repeat time for this orbit is 5.5781. It includes four symmetric turning points : 

{ q,P,C } = { 0.0000, ±1.5499, 0.0000 } and { ±1.2144,0.0000,0.0000 } . 

This orbit is the central core of an infinite family of nested tori. See Figure 6 . The tori are 
quasiperiodic regular structures tracing out two-dimensional private regions where there is 
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FIG. 5: Simplest [ shortest period ] periodic orbit for the Nose-Hoover oscillator with T = 1 . 
The initial condition is the turning point (q,p,() = (0.0,1.5499,0.0) . The motion equations are 
q = +p ; p = — q — (p ; C = P 2 ~ 1 • The ( q,p/s,(/ ) trajectory from the original Nose equations 
with s initially equal to unity and with # reduced from 2 to 1 [so that the motion equations are 
q = (p/s 2 ) ; p = —q ; C = (p 2 /s 3 ) — (1/s) ] traces out exactly the same ( q,p , £) trajectory pictured 
here, but with a period 7.1973 rather than 5.5781 . The Nose-Hoover Lyapunov exponent varies 
in the range ±0.6513 on this orbit. 

none of the chaos present which would make new three-dimensional regions accessible. Let 
us have a more detailed look at the mechanism enabling chaos by measuring the rates at 
which chaotic orbits separate on their exploratory journeys. 
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FIG. 6: Cross sections with q = 0, p = 0, ( = 0, starting with the periodic orbit penetration at 
(±1.2144,0,0) and increasing or decreasing the initial coordinate value of ±1.2144 by 15 successive 
steps of ±0.1, ending up with the initial q values ±2.7144 . The resulting cross-sections are shown. 
The 12th positive and negative initial conditions lie within the chaotic sea. All the rest generate 
tori which trace out individual points along the lines shown in this figure. 

C. Characterizing Chaos through Lyapunov Instability 

Chaos is an important topic. Without it there would be no hope for correspondence 
between the microscopic and macroscopic descriptions of material behavior. Chaotic tra¬ 
jectories exhibit “Lyapunov instability” — two nearby trajectories [ a “reference” and its 
“satellite”, as explained below ] rapidly separate. When averaged over the chaotic sea 
( that both of them inhabit ) the mean value of this instantaneous separation rate de¬ 
fines the largest ( positive ) Lyapunov exponent, Ai . This instantaneous separation rate , 
Ai(f) = | ri — r 2 |/| ri — r 2 | , is an observable which is easy to measure. 
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Giancarlo Benettin’s Lyapunov-exponent algorithm 1 - follows the simultaneous dynamics 
of a “reference” trajectory and a nearby “satellite”, rescaling their separation to maintain 
their closeness at the end of every timestep. Typically only the position of the satellite is 
adjusted to restore the length of the separation vector | r s — ry \ . Alternatively both of 
two neighboring trajectories rq and r 2 can be adjusted symmetrically. Evidently a theoreti¬ 
cal treatment becomes complicated at the boundary where toroidal and chaotic trajectories 
meet, mix, and coexist. A more elegant continuous version of Benettin’s step-by-step rescal¬ 
ing can be implemented by including a constraining Lagrange multiplier in the differential 
equations of motion^. In the symmetric case the multiplier is applied to both trajectories : 

h = v i + (A/ 2 )(r 2 - n) ; r 2 = v 2 + (A/2)(ri - r 2 ) . 

Here v represents the unconstrained motion equations while r — {q,p, C ) describes the actual 
constrained motion. The Lagrange multiplier enforcing the constraint of fixed separation is 
just the local value of the largest Lyapunov exponent. Its longtime average is the largest of 
the three global Lyapunov exponents which together constitute the “Lyapunov spectrum” 
of the three-equation Nose-Hoover oscillator : 

Ai = ( Ai (t) ) ; Ai (t) = (vi - v 2 ) • (ri - r 2 )/(ri - r 2 ) 2 . 

A 2 and A 3 are defined in terms of the comoving growth or decay rates of an area defined by 
three trajectories — the rate is Ai + A 2 ; and of an infinitesimal volume <g> defined by four 
trajectories with rate : 

(■ dq/dq ) + {dp/dp) + {d(/d() = 0 - ((t) + 0 = (<g>/<8>) = Ai(t) + A 2 (t) + A 3 (t) . 

For the isothermal Nose-Hoover oscillator with all the parameters equal to unity that equi¬ 
librium oscillator, chaotic or not, has no dissipation. This is implied by its close relationship 
to the four oscillator equations of motion according to Nose’s Hamiltonian formulation : 

q = {p/s 2 ) ; P= -q\ s = C; C = (T 2 / s3 ) - (W/s) . [ Nose ] 


Whether ff is chosen equal to 1 or to 2, the four Lyapunov exponents describing Nose’s 
oscillator sum to zero as a consequence of Liouville’s Theorem : 

{dq/dq) + {dp/dp) + {ds/ds) + {df/df) = (<g)/< 8 >) = \i(t) + A 2 {t) + A 3 {t) + A 4 (t) = 0 . 


16 


The mean value of the friction coefficient, ( ( ) , vanishes in both cases [ Nose or Nose- 
Hoover ] because the motion equations are “conservative” whether or not the initial condition 
is chaotic. Accordingly, although the local Nose-Hoover spectrum sums to — ((t) its time- 
averaged spectrum sums to zero : 

( (d\n®/dt) ) = ( (dp/dp) ) = (-() = Ai + A 2 + A 3 = 0 . 

In view of Hamiltonian mechanics’ time-reversibility both spectra are also symmetric , 
with the first and last time-averaged exponents summing to zero. Within the chaotic sea 
the Lagrange-multiplier analyses of the equilibrium oscillators give : 

{ A } = { Ai, A 2 , A 3 , A 4 } = { +0.001925, 0 . 000 , 0 . 000 , -0.001925 } ; [ Nose ] 

{ A } = { Ai, A 2 , A 3 } = { +0.0139, 0 . 000 , -0.0139 } . [ Nose - Hoover ] 

In the quasiperiodic toroidal regions where chaos is absent all the time-averaged Lyapunov 
exponents are zero. Although the time averaging produces symmetric spectra the instanta¬ 
neous spectra need not be symmetric. Let us demonstrate that somewhat surprising lack of 
forward-backward symmetry next. 

D. The Tricky Time Reversibility of the Nose-Hoover Lyapunov Spectrum 

One might well expect that the “local” ( instantaneous ) oscillator Lyapunov spectra 
are time-reversible too. After all, both the trajectories used to define the largest Lyapunov 
exponent are reversible. Growing separation, forward in time, corresponds to diminishing 
separation when reversed, and vice versa. A simple way to test this idea of time-reversible 
Lyapunov exponents is to store a reference trajectory { q,p, C } going forward in time, with 
dt > 0 . Then, analyzing the instabilities forward and backward in time, as described by 
the tendency of an adjustable “satellite” trajectory to flee or approach a stored “reference” 
is a fruitful approach. 

Time reversibility can be imposed on a stored reference trajectory { q,p, C } in either of 
two ways. The stored data can be used “as is” for a reversed trajectory simply by changing 
the sign of the time increment dt . Then the stored data are obviously solutions of the three 
reversed motion equations : 

q = —p ; p = +q + (p ; ( = 1 — p 2 for dt < 0 . 
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Alternatively one can simply change the signs of the stored { p } and { C } [ and { A } if the 
Lagrange multiplier approach is used ] . Then the modified sequence of { q,p, ( } satisfies 
the original motion equations when it is played backward. These two ways of exploring 
reversibility are precisely equivalent in computation. Only the signs of some of the numbers 
are changed. The magnitudes forward and backward are identical. 

If one were to follow this reversibility procedure for both the reference and the satellite 
trajectories, then increasing separation going forward in time would correspond to decreasing 
separation in the reversed trajectory, and hence to a negative Lyapunov exponent rather than 
a positive one. Evidently in the Nose-Hoover oscillator case, with just three exponents, the 
largest exponent going forward would become the most negative going backward : 

A{ = ( Ai (t) )dt >o = — {X$(t) )dt <o • 

Although the Lagrange-multiplier equations using a stored reference trajectory with stored 
satellite trajectories are all of them time-reversible, with both v(t) and Xi(t) changing sign 
if dt < 0 , the reversed exponents are erroneous ! The correct approach, storing a reference 
trajectory and then generating a new set of reversed-time satellite trajectories to go with 
the reference, using Lagrange multipliers or with Benettin’s rescaling algorithm, gives an 
instantaneous spectrum which is unrelated to its forward-in-time twin. For us this seems 
surprising, even though it is fairly well known, and suggests interesting research directions, 
perhaps leading to a better understanding of which spectra reverse and which do not. The 
simplest explanation is probably best : the tendency of two trajectories to separate can 
depend only on their past history, and not the future. On a particular time-reversible 
trajectory, { q(ndt ) } , without knowing whether dt is positive or negative, the notions 
of past and future could be thought nebulous. The low-cost readily-available cure for this 
uncertainty is straightforward computation and observation, finding out what does happen. 

What actually does happen is surprising, and is illustrated in Figure 7 . The largest 
Lyapunov exponent going forward in time and the largest going backward typically sum 
nearly to zero, indicating that the vectors joining the reference and satellite trajectories are 
almost parallel in the two time directions. The stretching or shrinking observed for dt > 0 
is replaced by its opposite, shrinking or stretching for dt < 0 : 

+Af >0 (t) ~ —Af <0 (t) . 
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FIG. 7: The summed spectra, both forward and backward, and the negative value of ( all 
correspond to the same bold curve in the lower right panel. The left panels represent the first 
two Lyapunov exponents forward ( thick lines ) and backward ( thin lines ) in time. The summed 
exponents, forward and backward for the two methods of time reversal, satisfy similar sum rules : 
Ai + A 2 + A 3 = ^ • The coordinate q (thick line) and the momentum p (light line) versus time are 
shown in the upper right panel. Notice that a brief segment of time, just past 2520, during which 
the momentum is near zero corresponds to the only part of the plotted interval in which the first 
Lyapunov exponents have very different magnitudes in the forward and backward directions. The 
arrow points out the maximum in A{ to which A^ is unrelated. 

Of course this cannot be precisely true as the averaged values of both versions of Ai are 
positive. But the fluctuations in the exponent are two orders of magnitude larger than the 
relatively-small time-averaged value of ±0.0139 . The fluctuations are of the order of the 
oscillator frequency, to = 2'nv , rather than the much smaller instability rate. 

A trajectory fragment [ 50,000 timesteps forward and backward from near the center 
a much longer simulation, 0 < t < 5000 ] and specially selected to show the reversibility 
phenomena, is analyzed in Figure 7 . Apart from the single strong maximum in Ai (t) dt>0 
indicated by the arrow the sum of the two exponents is close to zero. There is no apparent 
correlation between the values of the second or the third Lyapunov exponents in the two 
time directions. Thus trajectory stability depends qualitatively on the direction of time, 
even for equilibrium systems. 
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In any case, the sums of all three exponents going forward and backward in time both 
must obey the instantaneous sum rule : 

Ai(f) + A 2 (t) + A 3 (f) = (dln®/dt) = =fC • 

The equilibrium oscillator considered here behaves like a dissipative system with an asym¬ 
metry between the forward and backward directions of time which can be traced to the 
varying comoving phase volume <8>(£). With even a little thermostating the resulting heat 
transfer is enough to break the pairing symmetry expected for Hamiltonian systems. In 
nonequilibrium steady states we will see that this same symmetry breaking is unrelenting 
and in fact prevents reversing the formally time-reversible trajectories by any means other 
than reusing stored trajectories. 


IV. OSCILLATOR ERGODICITY VIA GENERALIZED FRICTION 

We have seen that numerical solutions of the three Nose-Hoover oscillator equations 
( QjPjC ) are f ar from ergodic. The same is true for the four Nose equations ( q,p,s,( ) 
because the two trajectories are identical for the corresponding scaled variables : 

tNH = St N -> { <?,P, (,t }nH = { Qi (p/ S), C, t }v( scaled ) • 

Provided that the initial values of { q,p,( } correspond at t NH = t N = 0 with s(t — 0) = 1 
plots of ((q) are identical for the two sets of differential equations. This equivalence is a 
useful demonstration of Nose’s “scaling of time” with the time-scaling variable s . 

Though ergodicity is lacking, the friction coefficient equations for £ guarantee that the 
second velocity moments are equal to the specified value of the temperature 

( p 2 ) NH = ( (p/s) 2 ) n = T. 

By controlling another moment the velocity distribution should come more closely to resem¬ 
ble the equilibrium Maxwell-Boltzmann Gaussian. In fact two moments can be enough 1 ^. 
Consider the simultaneous control of p 2 and p 4 . For this generalized Nose-Hoover oscillator 
problem two friction coefficients are involved, ( for p 2 and £ for p 4 : 

q = +p ; P = -q - Cp - £(p 3 /T) ; C oc (p 2 /T) - 1 ; £ oc (p 4 /T 2 ) - 3(p 2 /T) . [ HH ] 
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Solutions of these “Hoover-Holian” [ HH ] equations, with all parameters equal to unity, 
appear to be ergodic, with the probability density and all of its moments converging to 
Gibbs’ Gaussian moments better and better as the sampling time is increased. At the 
same time careful examination of the two-moment flow’s cross sections reveals that there 
are no “holes” in the doubly-thermostated chaotic sea. This is in marked contrast to the 
holy-sea sections seen in the singly-thermostated flows illustrated in Figures 4 and 6 . 
Comprehensive tests for ergodicity were formulated and applied to the doubly-thermostated 
[ HH ] equations in connection with the 2014 Snook Prized. 

Evidently one can never achieve perfect agreement in numerical tests as the probabilities 
of high-energy states are not only infinitesimal but would also require infinitesimal timesteps 
for computational stability. The four [ HH ] equations above, which control the second and 
fourth velocity moments, are only one of several methods for achieving ergodicity for the 
oscillator with the use of two friction coefficients. In view of this success it is natural 
to wonder whether or not a single carefully-chosen thermostating variable could make the 
dynamics ergodic. 

There are many different approaches to the ergodicity problem. The symmetry of the 
oscillator’s coordinate and momentum suggests that one could thermostat q 2 just as well 
as p 2 or perhaps even both. By interpreting q 2 as the fluctuation of the force these ideas 
can be, and have been, applied to more complicated systems. Extensions of this idea to the 
“weak” (time-averaged ) control of two or more different moments, with forces proportional 
to 


[ (f/T) - 1 ], [ (, J /T 2 ) - 3(, 2 /T) ], [ (p 2 /T) - 1 ], [ (p7T 2 ) - 3(p 2 /T) ],..., 

have proved fruitful. Rather than describe all of these efforts let us concentrate on the 
simplest of them, the simultaneous weak control of both ( p 2 /T) and ( p 4 /T 2 ) using a single 
friction coefficient. A useful tool in investigations of this sort is the y 2 test, which makes 
it possible to estimate numerically which of two different thermostat choices gives results 
closer to the Maxwell-Boltzmann velocity distribution. We turn to that next. 
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A. Ergodicity with a Single Thermostat Variable — the y 2 Test 


After 30 years of failed attempts involving dozens of investigators we recently discov¬ 
ered that ergodic control of a harmonic oscillator appears to be possible with just a single 
thermostat variable. An important clue, from the investigations of Bulgac, Kusnezov, Ju, 
and Bauer 21,22 , was that cubic thermostats provide enhanced ergodicity relative to linear 
ones. With cubic thermostats ergodic motion can occur even in the minimal case of a 
three-dimensional phase space. For maximum simplicity we assign any and all control vari¬ 
ables to the momentum because the canonical momentum distribution is always the same, 
independent of the chosen potential energy. 

Accordingly, let us assign the entire thermostating burden to the momentum through 
just two velocity moments, again using the harmonic oscillator to illustrate : 

q = p ; v = -q - C n [ ap + P{p 3 /t) + -y(p 5 /t 2 ) ] ; 

C = «[ o P 2 /T ) - 1 ] + /3[ 0 P 4 /T 2 ) - 3 (p 2 /T) } + 7 [ (p 6 /T 3 ) - 5 (p 4 /T 2 ) ] . 

With this functional form of control, where n = 1 or 3 , it is easy to show that the corre¬ 
sponding solution of Liouvillc’s continuity equation is again the Maxwell-Boltzmann Gaus¬ 
sian distribution : 

f(q,p,C) = e -« 2 / 2 T e-^/ 2 T e-^ +1 /(" +1 ) . 

Monte Carlo or grid-based exploration of (a, /3, 7 ) parameter space reveals many binary com¬ 
binations of two types, (a,/3, 0) and (a,0, 7 ) , which fill out Gibbs’ distribution. Candidate 
models have good velocity and coordinate moments as well as phase-space cross sections 
without discernable “holes”. Figure 8 shows ( = 0 cross sections for two successful combi¬ 
nations (a, (3, 7 ) = ( 0.05,0.32,0.00 ) and ( 1.50,0.00,-0.50 ). Because motion equations 
involving the fifth and sixth powers of velocity are “stiff” we have mostly restricted our 
detailed explorations to combinations of the second-moment and fourth-moment controls. 

How best to find such combinations ? There are many ways. A purely visual approach is 
relatively effective. A hundred-frame movie is a convenient visualization tool. The frames 
result from choosing an arbitrary grid of 100 ( a,/3 ) values and plotting 200,000 successive 
( q,p ) cross-section points for which ( vanishes. Frames lacking apparent holes in the 
density, and also providing good second, fourth, and sixth moments identify ( a, f3 ) choices 
as good candidates for ergodicity. The most interesting sections can be confirmed ergodic 
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FIG. 8 : Cross sections for ( a, /3 ,7 ) = ( +0.05, +0.32,0.00 ) [ at the left ] and ( 1.50,0.00, —0.50 ) 
[ at the right ]. These sections are for equilibrium oscillators [ T = 1 ] and use controls linear and 
cubic in p while linear in £ , with n = 1 . The white lines correspond to “nullclines” where the 
velocity normal to the cross section vanishes. The scales all range from —5.0 to + 5.0 . The upper 
panels are the £ = 0.0 sections. The lower panels are the q = 0.0 sections. 

with greater accuracy by zooming in, while using a few million section points rather than 

200,000 . 

A more systematic approach, also useful, but not at all foolproof, can be based on Pear¬ 
son’s x 2 statistic, x 2 comes in handy when it is desired to know if an observed distribution 
{ o } ( from a numerical simulation ) is consistent with a predicted one ( with expected values 
{ e } from a theoretical analysis of the flow equations ). For a candidate ( a, f3 ) combination 
a coarse-grained probability can be defined and determined within N discrete sampling bins. 
The mean-squared deviation of the bin probabilities from the expected theoretical value is 
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divided by the expected bin population and averaged over bins. The Central Limit Theorem 
predicts that the resulting sum : 

N 

X 2 = E( (°- e )V e ) , 

where the angle-bracket average is over the points in each of the bins, will approach the 
number of bins N if the “expected” already-known distribution { e } , is a match to the 
“observed” one, { o } . 

The simplest special case of this idea results if a uniformly flat distribution of random 
numbers { 0 < 72. < 1 } is divided into N equal sampling bins. With just two bins the 
large-sample limit is y 2 = 1 ; with four bins 3 ; with six bins 5, and so on. The limiting 
value of x 2 with N bins is just N — 1 . With data gathered from dynamical simulations 
where the distribution is not flat the dependence on the sample size is irregular and the 
convergence is slowed due to the inevitable serial correlation of sampled trajectory data. 

Knowing that Gibbs’ canonical oscillator distribution is Gaussian in all the variables 
makes it possible to test bins in 5 or p or (, or combinations of these variables, using the 
X 2 goodness-of-kt criterion. With 100 bins and a billion data points values of y 2 within 
ten percent of the number of bins, y 2 ~ 100 ± 10 , are typical when the distribution being 
observed really is Gaussian. It is possible to debug such a program using a good random 
number generator such as Press’ “ran2” generator from Numerical Recipes. 

Our investigations suggested that control variables based on the differential equations : 

q = p ■ p — —q — £”[ ap + P(p 3 /T) + 7 (p 5 /T 2 ) } ; 

C = «[ 0 P 2 /T) - 1 ] + (3[ 0 P 4 /T 2 ) - 3 (p 2 /T) } + 7 [ (p 6 /T 3 ) - 5(;p 4 /T 2 ) ] , 

can provide ergodicity for the harmonic oscillator. The equations are particularly stiff when 
the sixth moment is included so that (a,/3) = (0.05,0.32) is the most promising of the 
combinations we have tried to date. We have also included sample sections using a two- 
parameter (a, 7 ) model controlling the second and sixth velocity moments. Because these 
equations are stiffer, requiring timesteps of order 0.0001 rather than 0 . 001 , the second and 
fourth-moment control is the better choice. With this model the fluctuations of the largest 
local Lyapunov exponent correspond to a standard deviation of about 27, two orders of 
magnitude larger than its long-time-averaged value so that the exploration of phase space 
is relatively rapid. 
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Because the equilibrium Maxwell-Boltzmann momentum distribution e~ p2 ^ 2T is the same 
for any potential, this finding suggests that it may well be possible to thermostat any small 
Hamiltonian system in this way. Let us check this idea for the simple pendulum. 

B. Thermostating the Nose-Hoover Pendulum Problem 

We have seen that the Nose-Hoover oscillator can be thermostated in a wide variety of 
ways using either one, two ( or possibly even three ! ) control variables, though the stiffness 
suggests that using three controls is unwise. In each case the consistency of the solutions 
can be checked using Liouville’s Theorem to confirm that the stationary flow leaves Gibbs’ 
probability density unchanged : 


(df/dt) = -V • (fv) = 0 . 

To apply similar ideas to the pendulum problem^ using a single friction coefficient we 
need only to replace the potential energy: (q 2 / 2 ) —* — cos(g) so that Gibbs’ canonical 
distribution becomes : 

U = - cos(q) + (p 2 / 2 ) f(q,p, () oc e cos{q)/T e~ p2/2T e~ c+1/{n+l) with - n < q <+n . 

We can make the friction coefficient ( consistent with the Gaussian momentum distribution 
by using an arbitrary collection of moments, for instance : 

q = p ; p = - sm(g) - C n [ op + P(p 3 / T ) + 7 (p 5 /t 2 ) ] ; 

C = «[ (P 2 /T) -l]+/3[ 0P 4 /T 2 ) - 3 (p 2 /T) } + 7 [ (p 6 /T 3 ) - 5 (p 4 /T 2 ) ] . 

The cross sections for (a, (3) = (0.300, 0.300) [ found visually ] and (0.088, 0.188) [ found using 
Pearson’s y 2 test ] are shown in the Figure 9 and, at least from the visual standpoint, both 
distributions appear to be ergodic. The color indicates the magnitude of the local Lyapunov 
exponent, Ai(i) . 

Our oscillator and pendulum examples both suggest that there is a dynamics, determinis¬ 
tic and time-reversible, which closely follows Gibbs’ canonical distribution. Because many of 
the problems addressed with molecular dynamics involve isothermal rather than isoenergetic 
processes this makes isothermal molecular dynamics a particularly useful tool. 
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With weak polynomial control of the momentum it appears that we have completed Nose’s 
search for a deterministic time-reversible dynamics satisfying Gibbs’ canonical distribution. 
With equilibrium thermostats under control let us go on to consider the extension of this 
dynamics to nonequilibrium systems in which the temperature varies with location-. We 
demonstrate the possibilities by including a temperature gradient in the harmonic oscillator 
problem, allowing the dynamics to become dissipative by transferring heat from “hot” to 
“cold”. These problems have an intrinsic pedagogical interest because they are simultane¬ 
ously time-reversible and dissipative. They generate multifractal attractor-repellor pairs, 
often with a considerable aesthetic interest. 

V. NONEQUILIBRIUM TIME-REVERSIBLE DISSIPATIVE OSCILLATORS 

Gibbs’ equilibrium canonical distribution depends upon the temperature T — ( p 2 /mk ) 
where p is a cartesian momentum component for a particle with mass m . For simplicity we 
continue to choose both k and m equal to unity. We can introduce a nonequilibrium tem¬ 
perature gradient, VT, by choosing a coordinate-dependent temperature T(q) . This opens 
up the possibility for heat transfer leading to a quantitative treatment of nonequiibrium 
problems. We choose a smooth profile with a maximum temperature gradient e : 

T(—oo) = 1 — e < T(q) = 1 + etanh(g) < 1 + e = T(+oo) . 

Although the equations of motion , 

q = p ■ p=- q - (p- ( = C (p 2 ,p 4 ,p 6 ,T) , 

are still time-reversible [ with (+£, +q, +p, +£) —* (—t, +q, —p, —() } the dynamics can turn 
out to be dissipative and irreversible . How can this be? What does it mean? To address 
these questions, which are good ones, we must look at time reversibility in more detail. 

The concept of time reversibility^ 4 can be made unnecessarily complex by introducing the 
concept of phase-space involutions. A straightforward definition is the wiser choice: First, 
imagine a movie of the motion in question ( this presupposes a connection between the 
dynamical system of differential equations and objects capable of visual representation ) ; 
Second, play the movie backwards, ( but with the clock on the wall still recording a steady 
increase in the time ) just reversing the order of the frames. In the backward movie velocities 
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FIG. 9: Ergodicity of the isothermal Simple Pendulum with T = 1 and cubic control [n = 3 and 
p oc —C 3 ] • These cross sections correspond to (a,/3) = (0.300,0.300) to the left and (0.088,0.188) 
to the right. 7 = 0 in both cases.The vertical scales range from —4.0 to +4.0 and the horizontal 
scales range from —5.0 to +5.0 . Note that the range of the pendulum coordinate ( where q is an 
angle ) is periodic : — it < q < +n . 

change sign but coordinates do not. If the backward movie obeys the same equations as the 
forward one the dynamical system describing the motion is time-reversible. If not, then not. 
Variables odd in the time, such as velocity and the microscopic heat flux, change sign in the 
reversed motion, but parameters, like gravity, are held fixed. In every case that we study 
here our microscopic differential equations of motion satisfy this criterion. I 11 the simplest 
example, the Nose-Hoover oscillator , 

q = +p ; p = -q - (p ; C oc [ ( p 2 /T ) - 1 ] ; 
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reversing the sign of p and the variable £ is equivalent to reversing the sign of the time so 
that this system is time-reversible, even in the case where the temperature T depends on 
the coordinate q . In the simplest nonexample, 

q = +p ■ p = -q - (p , 

where ( is now a fixed parameter, not a variable, changing the sign of p corresponds to 
reversing the q equation but not the p one so that the constant-friction system is not time- 
reversible. The same is true of continuum solutions of viscous fluid flows and Fourier heat 
flow. The continuum constitutive relations for the shear stress of a Newtonian fluid with 
a velocity gradient and the heat current in a Fourier heat conductor with a temperature 
gradient : 

cr xy = v[ ipu x /dy) + ( du y /dx) ] ; Q x = -n(dT/dx) , 

are specially interesting. Here rj is the shear viscosity and k is the thermal conductivity. 
From the virial theorem we know that shear stress a is an even ( time-reversible ) function 
of the velocities, while the heat flux vector Q is odd. Both these observations contradict the 
phenomenological macroscopic constitutive relations laid down by Newton and Fourier. 

From thermodynamics we are well aware that “The entropy of the Universe increases.” 
Entropy is associated with heat reservoirs. When a reservoir absorbs heat SQ its entropy 
increases by 5Q/T . Likewise, when it releases heat the reservoir entropy decreases. Let us 
get back to the main question, “How do time-reversible motion equations produce irreversible 
behavior ?” 

It is a curious and hard-to-grasp fact that reversible mathematical equations can lead to 
irreversible behavior in the presence of Lyapunov instability, when the separation between 
two nearby trajectories increases. There would seem to be no reason why an increasing 
separation would win out over a decreasing one, particularly in the case where phase volume 
is conserved. But a positive Lyapunov exponent signals the system’s seeking out the direction 
of increased phase-space states. Think again of computing the largest Lyapunov exponent 
by rescaling the separation of two nearby trajectories. That is, consider the linear variation 
of a coordinate perturbation backward and forward in time, 5 , paralleling the direction 
characterized by the local Lyapunov exponent A (f) : 

5/8 = A(t) ( 5(t)/5( 0 ) ) = (1/2) [ e +Xdt + e~ Xdt } = cosh(A dt) , 



which gets bigger, indicating more phase volume, when averaged over the two possible 
time directions. In dynamical systems theory directions with growth are referred to as the 
“unstable” manifold while the decay directions are the “stable” one. If the direction of 6 
is allowed to develop “naturally”, with only its length, but not its direction constrained, it 
soon comes to point in the direction associated with the largest Lyapunov exponent Ai . 

But this is not at all the direction of maximum growth. To see this consider a simple 
example^, a harmonic oscillator’s orbit where the mass and the force constant are both equal 
to 1/4, and with an orbit perturbation ( 8q,8p ) having a fixed length | 5 | . If we choose 
a displacement parallel to the direction of the perturbation ( as in Benettin’s rescaling 
algorithm used to determine the maximum Lyapunov exponent ) we need to solve two 
coupled evolution equations : 

5 q = +4 5 p — X5g ; 8 P = —(1/4)5^ — A 8 P — > A = {15 / A)5 q 5 p /\ 5 | 2 , 

subject to the fixed-length constraint 8 q 5 q + 8 P 8 P = 0 . The maximum value, for equal 
perturbations in the two directions, gives a growth rate of (15/8) . On the other hand the 
unperturbed growth rate in a general direction, is (4 8 p , —(1/4 )8 q ) , which has its maximum 
value of 4 for a perturbation parallel to the p direction. Maximum and minimum growth 
rates are time-reversible but the dependence of the Lyapunov growth and decay rates are 
not. The latter rates depend upon past history, and not future destiny. 

A. Nonequilibrium Examples with Weak but Stiff Control of p 2 and p 6 

A specific ergodic system , 

q — p ; p = —q — ([ 1.5 p — 0.5 (p 5 /T 2 ) ] ; 

c = 1.5[ (p 2 /T) - 1 ] - 0.5[ (//T 3 ) - 5(p 4 /T 2 ) ] , 

exerts a passive control over the second and sixth velocity moments. Passive in that the 
equations are necessarily consistent with the Gaussian distribution but are not necessarily 
sufficient for that distribution to be realized. When the target temperature T(q) varies with 
coordinate , 

T(q) = 1 + etanh(g) ; 0 < e < 1 i —» 0 < T(q) < 2 , 
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the problem can become dissipative, with the comoving phase-volume rate-of-change nega¬ 
tive, 

( (®/<8>) ) = (dq/dq) + {dp/dp) + (<9(/<9() = ( -1.5C + 2.5( (p A /T 2 ) ) < 0 . 

Sample cross-sections of the corresponding chaotic sea appear in Figure 10. At equilibrium, 
not shown in the Figure, a ( q,p , 0) plot would show inversion symmetry, equivalent to noting 
that viewing the thermostated oscillator in a mirror would change the signs of (q,p) without 
changing the friction coefficient (. Changing the signs of the friction coefficient and both 
parameters (±£, ±1.5, ±0.5) leaves the equations of motion unchanged. 

An analytic calculation of the largest Lyapunov exponent for the (a, 7 ) problems illus¬ 
trated in Figure 10 requires the solution of the three coupled linearized equations of motion 
for the reference-to-satellite vector. In the p = 0 plane the equations for S — ( S q , S p , S ^ ) 
simplify : 

{ tig = d p - Ai Sq ; 6 p = -Sq - 1.5C<5 p - X 1 8 p ; S ( = -A^ c } . 

The constraint of constant length imposed by Ai in its role as a Lagrange multiplier, when 
applied in the p — 0 plane, makes it possible to relate the sign of Ai to that of the friction 
coefficient ( : 

(d/dt) [ 8^ + 8^ + S^ ] = 0 —> Ai = — a(d p ■ 

With a = 0.05 or 1.5 the upper half plane of the lower panels of Figures 10 and 11, 
with ( > 0 is entirely green, signifying a negative Lyapunov exponent. Positive values of Ai 
indicated by red in the Figures, correspond to negative values of Q in the p — 0 plane. As 
would be expected for a “friction” coefficient negative values promote growth and positive 
ones decay. As the temperature gradient increases the dissipation grows. Because the time- 
averaged dissipation is necessarily positive a comoving volume element d,qdpd( vanishes with 
time, exponentially fast. 

B. Nonequilibrium Example with Weak Control of p 2 and p 4 

Because the equations of motion for control of the sixth moment are stiff, requiring a 
Runge-Kutta timestep of order 0.0001 or 0.0002, we consider the simpler case of passive 
control of the second and fourth velocity moments : 

q — p ; p = —q — ([ 0.05p ± 0.32 (p 3 /T) ] 
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FIG. 10: Cross Sections of the Dissipative Oscillator with (a,j) = (1.5,—0.5) . For all of the 
plots the abscissa ranges from —5.0 to + 5.0 while the ordinate ranges from —4.0 to + 4.0 . 
The penetrations of the £ = 0 plane ( top panels ) corresponding to a positive largest Lyapunov 
exponent are in red, with negative values in green. The penetrations of the p = 0 plane ( shown 
in the lower panels ) are colored in the same way and confirm that negative values of the friction 
coefficient correspond to phase-volume growth. 

C = +0.05[ (p 2 /T) - 1 ] + 0.32[ (p 4 /T 2 ) - 3(p 2 /T) ] . 

Like those controlling ( p 2 ,p 6 ) these equations controlling ( p 2 ,p 4 ) are ergodic at equilib¬ 
rium, and are consistent with Gibbs’ distribution f(q,p, £) oc e~ q2 / 2T e~ p2 ! 2T e~^ 2 1 2 . 

Away from equilibrium, with T = 1 + etanh(g) and where tori are absent, these same 
equations generate dissipative fractal attractors. For small temperature gradients the di¬ 
mensions of the cross sections in Figure 11 and the dissipation vary smoothly while the 
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FIG. 11: Cross Sections of the Dissipative Oscillator with ( a, (5 ) = (0.05,0.32) and a friction 
linear in £ (n = 1) . For all of the plots the abscissa ranges from —5.0 to + 5.0 while the ordinate 
ranges from —4.0 to + 4.0 . 

motion remains ergodic. Just as in the Galton Board examples the resulting fractals signal 
the rarity of nonequilibrium states. 

Because the two-parameter motion equations are time-reversible there exists a symmetric 
set of ( = 0 cross-section states with reversed momenta ±p <—> =F p , reversed heat current, 
reversed dissipation, as well as phase-volume growth rather than collapse, all of these char¬ 
acteristics violating the Second Law of Thermodynamics. But all those states are unstable, 
with a reversed ( when time-averaged ) Lyapunov spectrum having a positive sum. They 
make up the unobservable and illegal repellor states. The repellor acts as a source and the 
attractor a sink for these time-reversible heat-flow problems. Despite their symmetry the 
repellor and attractor measures are different, zero and one, respectively. 
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FIG. 12: The effect of dissipation on phase-space cross sections for the chaotic, but noner- 

godic, Nose-Hoover oscillator. Sections with vanishing momentum are shown for four values of e 
: 0.00, 0.10, 0.20, 0.30 . The equations of motion are q = p ] p = — q — (p ', ( = [ ( p 2 /T ) — 1 ] , 
where T = 1 + e tanh(g) . 

The implication of Figure 11 is relatively simple. Starting out with an equilibrium dy¬ 
namics which follows Gibbs’ canonical distribution ( without holes ) relatively simple mul¬ 
tifractal attractors respond to a thermal gradient. Increasing the temperature gradient 
leads to a reduced attractor dimensionality and increased dissipation. The general approach 
of perturbing an ergodic equilibrium Gibbs’ ensemble evidently leads to relatively simple 
nonequilibrium steady states. In an effort to see whether or not this simplicity has a coun¬ 
terpart in nonergodic Hamiltonian mechanics we return to the Nose-Hoover oscillator and 
expose it to a thermal gradient next. 
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C. Nonequilibrium Nose-Hoover Oscillator Dynamics 


Although the equilibrium Nose-Hoover oscillator has all the complexity of Hamiltonian 
mechanics - chains of islands and their elaborations, it is worth exploring whether or not 
any simplification results away from equilibrium. Our previous work” indicated a crude 
boundary between small-gradient and large-gradient behavior around e = 0.40 . In Figure 
12 we consider strange-attractor solutions corresponding to three values of e, 0.10, 0.20, 0.30 
in addition to the equilibrium case. We use nonequilibrium versions of the Nose-Hoover 
oscillator with ( = [ (p 2 /T) — 1 ] rather than C = [ P 2 ~ T ] because the multimoment 
models are much simpler to formulate when the distribution of the friction coefficient (s) is 
not explicitly temperature dependent. 

We considered relatively long runs ( 10 12 timesteps ) in an effort to assure convergence. 
A bit more than halfway through the e = 0.10 simulation, with a fourth-order Runge-Kutta 
timestep of 0 . 01 , the chaotic strange attractor suddenly began to generate a torus which 
then gradually shrunk with time. Was this real, or an artefact? A check calculation with 
fifth-order Runge-Kutta, also using a timestep of 0.01 revealed no such behavior, exhibiting 
instead a chaotic solution. This problem illustrates the virtue of comparing results from the 
two or more integrators, particularly when longer runs are desirable. A careful investigation 
shows that the single-step fourth-order and fifth-order errors: 

RK4 error ~ —dt 5 /120 ; RK5 error ~ +dt 6 /720 , 

are in opposite directions, with the fifth-order still noticeable with a timestep dt = 0.01 in 
double-precision calculations. With dt = 0.001 both the fourth-order and fifth-order errors 
are negligible in double-precision work. For the reader interested in exploring these small 
effects^ an initial condition very close to the border between chaos and tori is (q,p,() — 
(<5, S , 3) with 6 small. A small nonzero value of 5 ( such as 10 ~ 12 ) is necessary to avoid the 
analog of a ( q,p ) fixed point in the (q,p,( = 0 ) cross section. 

Figure 12 shows the sign of the local Lyapunov exponent Ai(t) in color, both at and away 
from equilibrium ( e = 0.0 to 0.3 ). Notice that the near inversion symmetry in the (q,p, 0) 
plane for e = 0.10 , gives way to predominating Lyapunov instability far from equilibrium, 
at e = 0.30 . Below, the (g, 0, () plane shows that the sign of the local Lyapunov exponent 
is a perfect match of the sign of the local friction coefficient. Simply reversing the direction 
of the flow in the three nonequilibrium panels, corresponding to reflection of p about the q 
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axis, might be expected to change the signs of the Lyapunov exponents, but even close to 
equilibrium this does not happen. This is because that exponent depends upon the past so 
that there is a fundamental lack of symmetry in the local exponents. 

VI. LIOUVILLE’S THEOREM APPLIED TO NONEQUILIBIUM FLOWS 

We have seen that the continuity equation is an invaluable tool in finding constrained 
dynamical systems consistent with Gibbs’ canonical ensemble. This idea was used by Green 
and Kubo to express transport coefficients in terms of equilibrium fluctuations. Nonequilib¬ 
rium simulations, even far from equilibrium, use this same tool. Because the Galton Board 
as well as all of the thermostated oscillator problems we have considered involve three- 
dimensional flows it is natural to consider their analysis and display from the standpoint of 
Liouvillc’s phase-space continuity equation : 

(df/dt) = -V • (fv) i—¥ f = ( df/dt ) + v • V/ = -/ • Wv , 

also in three dimensions. The corresponding motion equations are represented by v = 
( q,p,C ) in the oscillator problems. 

It is tempting to imagine solving the flow equation directly, replacing the derivatives 
(q,p, C) by finite differences. Our colleague John Ramshaw made us the welcome present of 
his “upwind-differencing” computer program, which transfers density across all six faces of 
each cubic computational cell according to the velocities evaluated at the cell boundaries, 
determining which of any two adjacent cells is the “donor” of probability, and which is the 
“acceptor”. The flow is taken to be proportional to the donor probability though it would 
appear that an average probability is nearly as plausible. Using the average, however, leads 
to exponential instability 27 . At equilibrium, or in a steady state, the six flows into and 
out of every cell must balance. Evidently the algorithm conserves probability exactly, but 
not time-reversibly. Liouville’s continuity equation is time-reversible. But it is easy to see 
that in a reversed implementation of the flow algorithm the cells furnishing the probability 
forward in time will not have it returned exactly in the “reversed” step. 

The time-reversal thought experiment ( as well as its computational realization ) show 
that the algorithm has a diffusive irreversibility. To see this begin with a single filled zone. 
Going forward in time this zone would typically donate probability to half (three) of its 
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neighbor zones. The other three neighbors would remain empty. On reversal most of the 
probability gained by the three neighbors would remain there. Such a forward-plus-backward 
two-timestep simulation reveals why it is that a straightforward application of Liouvillc’s 
flow equation fails. 

A simple computational test of this Liouville algorithm is to reproduce the periodic 
conservative Nose-Hoover orbit, shown in Figure 5 . The initial conditions (q,p, () = 
(0.00,1.55, 0.00) generate it. Solving the Liouville equation with all of the initial density 
in the cells nearest the initial condition leads to an amplitude loss proportional to the cell 
size. After a few such damped oscillations the final quiescent stationary state is obtained, 
with most of the density bunched near the origin. Figure 13 displays the mean values of 
the finite-difference algorithm’s history for ten oscillation periods. Limit-cycle and chaotic 
problems lead to very similar results. 

A better rendition of the continuum motion equations could follow a set of Lagrangian 
points, attracted toward one another with a Lagrange multiplier designed to follow the co¬ 
moving flow volume precisely. Similar algorithms, based on interface tracking could perhaps 
be developed, but in the end the numerical solution of Liouville’s continuity equation appears 
much better suited to developing consistent motion equations than to evolving a continuous 
phase-space density. 

VII. ANALOGIES WITH MANY-BODY PROBLEMS 

Over the years thousands of papers have described the use of time-reversible control 
variables to solve simple problems. Shear flows, heat flows, and shockwaves are examples 
that spring to our minds. Large-scale biomolecular simulations use exactly these same 
ideas. In stationary nonequilibrium flows the atomistic forces are supplemented by boundary, 
constraint, and driving forces in such a way as to generate a nonequilibrium steady state. 
In every case the resulting phase-space distribution is fractal, representing a flow from a 
strange repellor, with a positive summed Lyapunov spectrum to a strange attractor, with 
a negative Lyapunov sum. Although the fundamentals are no more complicated than the 
examples detailed in this paper, the possibilities for more complex applications are and will 
continue unlimited. 
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FIG. 13: The Nose-Hoover periodic orbit with T = 1 , as in Figure 5 , is shown in red, with a 
bullet indicating the initial condition, (q,p,() = (0.00,1.55,0.00) . The finite-difference Liouville 
continuity equation solution for ( p(( q )) ) is shown for a time of 55.78, ten oscillation periods, in 
black, using a 40 x 40 x 40 mesh with a mesh spacing of 0.1 . Solutions with three finer meshes ( 
0.0500, 0.0250, 0.0125 ) establish that the artificial damping is proportional to the mesh size as 
is shown in John Ramshaw’s book, Elements of Computational Fluid Dynamics —. 
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